DeriveCurvature Subroutine

public subroutine DeriveCurvature(dem, curvature, profileCurvature, planCurvature)

Calculates the curvature [1/m] of a raster surface, optionally including profile and plan curvature. If any neighborhood cells are NoData, they are first assigned the value of the center cell. It calculates the second derivative value of the input surface on a cell-by-cell basis. For each cell, a fourth-order polynomial of the form:

 Z = Ax²y² + Bx²y + Cxy² + Dx² + Ey² + Fxy + Gx + Hy + I

is fit to a surface composed of a 3x3 window. 

 A = [(Z1 + Z3 + Z7 + Z9) / 4  - (Z2 + Z4 + Z6 + Z8) / 2 + Z5] / L4
 B = [(Z1 + Z3 - Z7 - Z9) /4 - (Z2 - Z8) /2] / L3
 C = [(-Z1 + Z3 - Z7 + Z9) /4 + (Z4 - Z6)] /2] / L3
 D = [(Z4 + Z6) /2 - Z5] / L2
 E = [(Z2 + Z8) /2 - Z5] / L2
 F = (-Z1 + Z3 + Z7 - Z9) / 4L2
 G = (-Z4 + Z6) / 2L
 H = (Z2 - Z8) / 2L
 I = Z5

   curvature = -2 (E + D)

   profileCurvature = -2 (D G^2 + E H^2 + F G H) / (G^2 + H^2)

   planCurvature  = -2 (D H^2 + E G^2 - F G H) / (G^2 + H^2)

Profile curvature is parallel to the direction of the maximum slope. A negative value indicates that the surface is upwardly convex at that cell. A positive profile indicates that the surface is upwardly concave at that cell. A value of zero indicates that the surface is linear. Planform curvature (commonly called plan curvature) is perpendicular to the direction of the maximum slope. A positive value indicates the surface is sidewardly convex at that cell. A negative plan indicates the surface is sidewardly concave at that cell. A value of zero indicates the surface is linear. Plan curvature relates to the convergence and divergence of flow across a surface.

References:

Moore, I. D., R. B. Grayson, and A. R. Landson. 1991. Digital Terrain Modelling: A Review of Hydrological, Geomorphological, and Biological Applications. Hydrological Processes 5: 3–30.

Zeverbergen, L. W., and C. R. Thorne. 1987. Quantitative Analysis of Land Surface Topography. Earth Surface Processes and Landforms 12: 47–56.

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(in) :: dem
type(grid_real), intent(out) :: curvature
type(grid_real), intent(out), optional :: profileCurvature
type(grid_real), intent(out), optional :: planCurvature

Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: D
real(kind=float), public :: E
real(kind=float), public :: F
real(kind=float), public :: G
real(kind=float), public :: H
real(kind=float), public :: L
integer, public :: i
integer, public :: j
real(kind=float), public :: z1
real(kind=float), public :: z2
real(kind=float), public :: z3
real(kind=float), public :: z4
real(kind=float), public :: z5
real(kind=float), public :: z6
real(kind=float), public :: z7
real(kind=float), public :: z8
real(kind=float), public :: z9

Source Code

SUBROUTINE DeriveCurvature &
!
(dem, curvature, profileCurvature, planCurvature)


IMPLICIT NONE

!arguments with intent in
TYPE(grid_real),    INTENT(in):: dem

!arguments with intent out
TYPE (grid_real), INTENT (out) :: curvature
TYPE (grid_real), OPTIONAL, INTENT (out) :: profileCurvature
TYPE (grid_real), OPTIONAL, INTENT (out) :: planCurvature

!local variables
INTEGER :: i,j
REAL (KIND = float) :: D, E, F, G, H
REAL (KIND = float) :: z1, z2, z3, z4, z5, z6, z7, z8, z9
REAL (KIND = float) :: L

!-----------------------------end of declarations------------------------------

!allocate grids
CALL NewGrid (curvature, dem)

IF (PRESENT(profileCurvature)) THEN
  CALL NewGrid (profileCurvature, dem)
END IF

IF (PRESENT(planCurvature)) THEN
  CALL NewGrid (planCurvature, dem)
END IF


!loop through dem grid
idim: DO i = 1, dem % idim
  jdim: DO j = 1, dem % jdim
    IF (dem % mat(i,j) /= dem % nodata) THEN
    
       !set length
       L = CellArea (dem,i,j) ** 0.5
       
       !set z5
       z5 = dem % mat (i,j)
    
       !set z2
        IF ( .NOT. IsOutOfGrid (i-1,j,dem) ) THEN
          IF ( dem % mat (i-1,j) /= dem % nodata) THEN
             z2 = dem % mat (i-1,j)
          ELSE
             z2 = dem % mat (i,j)
          END IF  
        ELSE 
           z2 = dem % mat (i,j)
        END IF
        
        !set z3
        IF ( .NOT. IsOutOfGrid (i-1,j+1,dem) ) THEN
          IF ( dem % mat (i-1,j+1) /= dem % nodata) THEN
             z3 = dem % mat (i-1,j+1)
          ELSE
             z3 = dem % mat (i,j)
          END IF  
        ELSE 
           z3 = dem % mat (i,j)
        END IF
        
        !set z6
        IF ( .NOT. IsOutOfGrid (i,j+1,dem) ) THEN
          IF ( dem % mat (i,j+1) /= dem % nodata) THEN
             z6 = dem % mat (i,j+1)
          ELSE
             z6 = dem % mat (i,j)
          END IF  
        ELSE 
           z6 = dem % mat (i,j)
        END IF
        
        !set z9
        IF ( .NOT. IsOutOfGrid (i+1,j+1,dem) ) THEN
          IF ( dem % mat (i+1,j+1) /= dem % nodata) THEN
             z9 = dem % mat (i+1,j+1)
          ELSE
             z9 = dem % mat (i,j)
          END IF  
        ELSE 
           z9 = dem % mat (i,j)
        END IF
        
        !set z8
        IF ( .NOT. IsOutOfGrid (i+1,j,dem) ) THEN
          IF ( dem % mat (i+1,j) /= dem % nodata) THEN
             z8 = dem % mat (i+1,j)
          ELSE
             z8 = dem % mat (i,j)
          END IF  
        ELSE 
           z8 = dem % mat (i,j)
        END IF
        
        !set z7
        IF ( .NOT. IsOutOfGrid (i+1,j-1,dem) ) THEN
          IF ( dem % mat (i+1,j-1) /= dem % nodata) THEN
             z7 = dem % mat (i+1,j-1)
          ELSE
             z7 = dem % mat (i,j)
          END IF  
        ELSE 
           z7 = dem % mat (i,j)
        END IF
        
        !set z4
        IF ( .NOT. IsOutOfGrid (i,j-1,dem) ) THEN
          IF ( dem % mat (i,j-1) /= dem % nodata) THEN
             z4 = dem % mat (i,j-1)
          ELSE
             z4 = dem % mat (i,j)
          END IF  
        ELSE 
           z4 = dem % mat (i,j)
        END IF
        
         !set z1
        IF ( .NOT. IsOutOfGrid (i-1,j-1,dem) ) THEN
          IF ( dem % mat (i-1,j-1) /= dem % nodata) THEN
             z1 = dem % mat (i-1,j-1)
          ELSE
             z1 = dem % mat (i,j)
          END IF  
        ELSE 
           z1 = dem % mat (i,j)
        END IF
        
        !compute D, E, F, G, H
         D = ((z4 + z6) /2. - z5) / L**2
         E = ((z2 + z8) /2. - z5) / L**2
         F = (-z1 + z3 + z7 - z9) / (4. * L**2)
         G = (-z4 + z6) / (2. * L)
         H = (z2 - z8) / (2. * L)
         
         curvature % mat (i,j) = -2. * (E + D)
         
         IF (PRESENT(profileCurvature)) THEN
           IF ( (G**2. + H**2.) > 0.) THEN
               profileCurvature % mat (i,j) = -2. * (D*G**2. + E*H**2. + F*G*H) / (G**2. + H**2.)
           ELSE
               profileCurvature % mat (i,j) = 0.
           END IF
         END IF
         
         IF (PRESENT(planCurvature)) THEN
            IF ( (G**2. + H**2.) > 0.) THEN
               planCurvature % mat (i,j) = -2. * (D*H**2. + E*G**2. - F*G*H) / (G**2. + H**2.)
            ELSE
                planCurvature % mat (i,j) = 0.
            END IF 
         END IF
       
    
    END IF
  END DO jdim
END DO idim


RETURN

END SUBROUTINE DeriveCurvature